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Abstract 


^We-descr-ibe a prototype network planning model for the U.S. Air 
”■ c. 

Traffic control system/ The model encompasses the dual objectives of 
managing collision risks and transportation costs where traffic flows 
can be related to these objectives. The underlying structure is a 
network graph with nonseparable convex costs; the model is solved effi- 
ciently by capitalizing on its intrinsic characteristics. Two specialized 
algorithms for solving the resulting problems are described: (1) truncated 
Newton, and (2) simplicial decomposition. 

The feasibility of the approach is demonstrated using data collected 
from a control center in the Midwest. Computational results with different 
computer systems are presented — including a vector supercomputer 
(CRAY-XMP) . The risk/cost model will have two primary uses: (1) , as a 
strategic planning tool using aggregate flight information, and (2) as an 
integrated operational system for forecasting congestion and monitoring 
(controlling) flow throughout the United States. In the latter case, 
access to a supercomputer will be required due to the model’s enormous 
size. 
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1. Introduction 

By the year 2000, the U.S. Federal Aviation Administration (FAA) 
plans to spend over $11 billion oh a program to upgrade and consolidate the 
U.S. air-traffic system — called the National Airspace System Plan [NASP]. 

On an average day in 1984 approximately 17,000 aircraft traveled these 
routes. It is is expected that traffic volume will increase by 62% over 
the next decade. The FAA has responsibility for managing the air-traffic 
system, especially those elements affecting risk. Indeed, the agency has 
been able to reduce overall (and relative) risks during the past, thirty 
years [ 1 ]. As new hardware and software technologies have been implemented, 
the FAA has provided mechanisms for monitoring factors associated with 
system risks. Recently a trend to deregulate the airline industry in the 
interest of greater efficiency and competitiveness has begun. Despite this 
trend the FAA must possess efficient procedures for assessing (and ultimately 
controlling) system risks. 

In this report we discuss the development of a prototype network 
planning model for flights on high-altitude jet routes over the U.S. air- 
space, in conjunction with the NASP. The model encompasses the dual objective 
of assessing risk-related measures and transportation costs. The underlying 
mathematical model has the special structure of networks graph. Since safety 
components depend upon interacting variables, the proposed model falls in 
the category of a nonlinear network with nonseparable cost functions. 

The model can serve as the building block for a management Information 
system that can assist FAA in two basic settings: 

(a) As a strategic planning tool, or 

(2) As an operational planning program for air-traffic 
scheduling and routing. 
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In the former case it is essential that the optimization can be solved 
efficiently and accurately. In this way extensive sensitivity analyses 
can be carried out and answers to "what if” questions obtained in a 
comparative fashion and without regard to numerical instability. The 
use of the model as an operational planning tool depends on the ability 
to "solve" the model under real-time conditions. 

By capitalizing on the special structure of the underlying network 
basis (a forest of 1-trees) we have developed highly efficient solution 
algorithms [2, 12, 13]. Thus large-scale problems can be handled within 
minutes or hours; repeated runs can be made and the use of the model for 
planning is feasible. Going a step further we have streamlined our 
algorithms for the architecture of vector supercomputers, in particular 
the CRAY-XMP at Boeing Computer Services. By taking advantage of the 
Cray's pipeline features we were able to solve the test cases in a matter 
of seconds — thus demonstrating feasibility of the network model as an 
operational planning tool. 

The rest of this paper is organized as follows. The mathematical 
models are defined in Section 2. A general description of two solution 
algorithms — ' (1) Truncated Newton, and (2) simplicial decomposition — 
appears in Section 3. Data requirements for the model are listed in 
Section A. Next, Section 5 presents the modeling of a representative 
air-traffic control sector (Indianapolis); some computational 
results are given in Section 6. Finally, we discuss limitations of the 
current model and directions for future research. 
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2. Air Traffic Control ModelinR 

A stochastic programming model for the air-traffic control problem 
was proposed in the early work by Ferguson and Dantzig [5]. More recently, 
general aspects of air-traffic systems planning and design were the focus 
of a conference at Princeton University [1 ]• We propose here a network 
formulation for the air-traffic planning problem. 

The generic nonlinear network model takes the form: 

[NLGN] Minimize F(x) 

Subject to: 

-ij - J.-Vkl ' ‘'l- f'"' 1 e N (2.1) 

j ei , ke6 . 

1 1 

1.. < X.. < u.., for (i,j) e E (2.2) 

ij - ij - ij ’ 


where 


F(x) = 
{N} = 
{E} = 


^ « t 

ij 

“ki 



convex function 
set of nodes 
set of arcs (edges) 
flow over arc (i,j) 
multiplier* on arc (k,i) 
supply/demand for node i 
I (i,j) c E) 

{ji(i.j) e E) 

{j|(j,i) e E) 


u..(f. .) upper (lower) bound on arc (i,j) 

X = {x..|x.. satisfies constraints (2.1) and (2.2)) 

13 ' 13 

*Multipliers are indicated so that airport congestions and other factors 
can be modeled. These features will be descussed In Section 7. 
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We may rewrite [NLGN] in the following form: 

[NLGN] Minimize F(x) 

Subject to X e X 
where 

X = {x]a • X = b.ii ^ u) 

Since the [NLGN] basis is a ’’forest of 1-trees" (collection of subtrees 
with one extra arc per subtree creating a single loop), efficient procedures 
are available for storing and updating the basis and other aspects of the 
algorithm [9, 12]. 

Figure 1 depicts a simplified example graph of the first planning model. 

This example consists of five airports and interconnecting routes, representing 

an aggregate of individual flights into and out of the designated airports. 

Note that each airport has three triangular nodes, indicating net traffic, 

number of incoming flights and the number of outgoing flights. Network arcs 

points in the obvious direction of traffic movements. 

The corresponding optimization model is shown below: 

[OPTl] Minimize {w, * [ I \ c^,z5.] + w„* s(z)} = F(*) 

^ (i,j)CA^ rCR^j ij ij ^ 

Subject to: 


X 


i 


I 

k£L 


ik 


iCA 


^i “ ^ ^ki 

kEM. 

1 




\ A 

"'ij - ""ij 

(i.k)eAj 


^kj - \j 

(k,j) such that jEA,k£M^ 

ik = 

I ^ik 

iEA, kEL^ 


livnt IfinttofiA 


5- 


ORiQiNAL PAQE 
OF POOR QUALITY 



Figure 1; Example Air-Traffic Network (Strateg ic Model) 
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0 < z^. < uT. 
- - iJ 


jEA.kCM^ 

reR^. and (i.j)eAj^ 



i£A 


W1.W2 ^ 0 


where we have used the following decision variables: 

x^(y^) = number of total outbound (inbound) flights, airport iCA 
= number of outbound flights from airport i to airport k 

z^j = number of flights over route r, airport i to airport j 

= number of inbound flights from airport k to airport i. 

Notation: 


A 

(i. j) 



set of airports 

feasible pair (linking airport i to airport j) 
all feasible routes for pair (ij) 


M. 


'ij 


b. 

1 


F(z) 


outbound destinations for airport i 
inbound airports for airport j 
{(i, j) I i£l and 

marginal transportation cost on route 
net traffic at airport i 

convex nonseparable function, measuring risk 
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This single period planning model aggregates all individual flights 
between airport-pairs over the planning horizon. While extending the model 
to a large number of multiple periods will considerably increase the model's 
size, it would be, conceptually, straightforward to depict several periods. 

As such, the [OPTl] model can be used for supporting airport resource 
decisions, such as opening new runways and expanding the capabilities of 
control towers. The primary aims are to identify bottlenecks and to predict 
imbalances for the U.S. air-traffic system, given a variety of scenarios. 

In this representation the optimization model takes a bi-criteria 
objective function consisting of risk and systemwide cost, with respective 
weights w^ and w^ . It should be stressed that the weights are used only as 
part of the sensitivity analyses to identify the efficient frontier and are 
not meant to be set a priori . Transportation cost is easily quantified in 
terms of traveled distance and fuel burn rates of the different aircraft 
models. A monetary value can also be assigned on factors like customer 
dissatisfaction due to flight delays and so on. Quantifying the risk 
component of the system, however, is a critical and controversial issue. 

At best we may consider optimizing a relative risk measure: the key idea 
is to compare some risk norm under different scenarios or system states, 
or with other similar systems. Odoni and Endoh [15] consider a probabilistic 
analysis of risk. In Section 5 we touch upon the issue of modeling risk as 
a function of congestion in the target sector during time intervals of 


interest . 



- 8 - 


The issue of identifying an acceptable point on the efficient risk- 
cost frontier is complex. It is obviously difficult for society to 
compare human risks and congestion delays against financial criteria. 

See, for example, Rowe [18]. Ultimately the solution depends upon societal 
tradeoffs and can only be derived through an informed political process. 

The network framework provides essential input data for this process by 
tracing out the efficient frontier; see Section 5. 

The aggregate model [OPTl] is unable to provide adequate details for 
operational planning purposes. 

Therefore, a second network model has been developed and an example 
is shown in figure 2. Here, decision variables monitor possible delays and 
alternative altitudes for every flight departing to or arriving at the 
airports of interest. The corresponding mathematical respresentation is 
defined below. In addition to optimizing transportation and delay costs 
this model limits the number of flights traversing regions of interest — 
e.g., control sectors — as a surrogate for miniiiizing risk. Including 


these aspects causes [0PT2] to be considerably larger than [OPTl]. 


[0PT2] Minimize T T T c?^ 
tir pcP^ teij. * 


fcF ttT, 





for all fEF 
for all feF, teT^ 


for all aCA, kEK 


o> 
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r a a a 

\ i \ i “k 


for all a£A, k£K 


y for all reS, keK 

^ I “ r 




where we have used the following decision variables: 


xf - 0.1 


dj = 0,1, tcT^ 


wf for keK 


y,, h ~ 1,2,...,K 

iC ni3x 


flight f follows route p during time period t 

flight f delays its departure during period t 

number of planes delayed at destination 
airport a during period k 

number of planes landing at destination 
airport a during period k. 


Notation: 


acA : set of airports 
feF : set of flights 




set of possible routes for flight f 


k = 1,2,...,K^^^ : time of periods in planning horizon 

tjETj : set of possible time periods for departure of flight f 

; set of flights traversing control sector r at period k 

: set of flights whose destination is airport a, 
and arrival time is k 

c^^ : marginal transportation cost for flight f on route 
p during period t 

c^ : marginal departure delay cost for flight f during 
period t 
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c^, : marginal landing delay cost at airport a 

during period k 

3d 

: upper (lower) limit on the total number of 
flights landing at airport a during period k 

It 

: limit on the number of flights traversing control 
sector r during period k 

S : of control sectors 

This model controls individual flights: for every flight it specifies 

possible altitudes, routes to be followed and any departure/landing delays. 

The objective function incorporates the marginal cost of transportation, 

and associated delay costs.* Transportation risk is assessed by imposing a 

limit on the number of flights through a particular sector during every time 

interval of interest. As mentioned this formulation provides more details 

than [OPTl] since it controls flights on an individual basis. The added 

details result in a more complex model: the model is not only larger than 

[OPTl] but also it has integer variables. It is, however, equivalent to 

[OPTlJin the risk minimization aspects where congestion in an identified 

region provides a surrogate for risk. By varying the limit on the number of 

flights traversing a control section during every time period we may achieve 

the same result as by varying the weights in the multi-objective formulation 

in [OPTl]. Mathematically the first formulation is more tractible: it deals 

with a continuous nonlinear network model, while the second model deals with 

a lihear multicommodity network problem with integer variables. 

*Note that the model minimizes total transportation costs, following a 
utilitarian economies. If this approach adversely impacts individual carriers, 
the model can be adjusted through multiple weightings or by means of side 
payments between carriers (or the FAA) . 
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The two models are expected to serve different purposes. The 
strategic planning model is intended to assess aggregate data, while the 
operational planning model deals with more microscopic aspects of air- 


traffic control. 



- 13 - 


3. Solution Algorithms 

We describe here the main features of two nonlinear programming 
algorithms that have been specialized to solve [NLGN] . The methods are: 

(1) A second order algorithm based on truncated Newton search directions 
(TN) , and (2) A first order algorithm based on the simplicial decomposi- 
tion of the feasible region (SD) . Both algorithms possess distinct 
characteristics: the truncated Nevrton algorithm can identify solutions 
to a high degree of accuracy; the simplicial decomposition algorithm can 
quickly converge to an approximate solution- The interested reader should 
refer to the papers by Ahlfeld et al. [2] and Mulvey et al. [13] for more 
detailed description of these algorithms and some computational results. 

3.1 Truncated Newton Algorithm 

In a manner similar to most nonlinear programming procedures, each 
(TN) iteration consists of two stages: (1) a search direction routine, and 

(2) a step length routine. Table 1 depicts the overall flow, in which the 
notation refers to the model [NLGN] presented in Section 2. 

The search direction must fulfill certain essential features so that 

the overall algorithm will converge and so that performance efficiencies 

are attained. First, the direction must both maintain feasibility and 

point downhill (in a minimizaxion context). Defining the search direction 

“k til 

as p and given a feasible point x at the k iteration, the usual Newton 

-k 

method for calculating p would solve the following quadratic programming 
problem: 
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Minimize <■' I, -k.. t„.-k,-k , ,-k^t-k 

[QP] ) G(x )p + g(x ) p 

P 


Subject to: 


X 

A • p 

p ^>0 

pJiO 


if X. 

3 

. , k 

if X . 

3 




= u , 


where 

g(x^) gradient of F(x^) at x^ 

G(x^) Hessian of F(x^) at 

By restricting our attention to a special projected matrix 2, whose columns 
form a basis for the null space of A, i.e., A • 2=0, the problem [QP] can 
be solved using the following two formulae: 


(2^ G 2). p^ = 2^^i^ 
-k „t -k 

P = Z pg 


(3.1. i) 


where 


2 = 


-b‘^s 

s 

N 


m 

s 

n-s-m 


-k 


G semi-positive approximation to the Hessian at point x 

-k -k 

g gradient of objective function at point x 

and where the decision variables have been partitioned into three sets: 

A = [b|s1n] 


g(x) = Igblgs'Sn’ 


-k 

P 


fPblPs'Pn^ • 
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The benefits of the Newton direction p are greatest in the neighbor- 
hood of a solution; however, it is expensive to calculate the solution of 
equation (3.1.i). In response, we adjust in a dynamic fashion the degree 
of accuracy of solving [QP]. A forcing sequence {tiM — 0 is employed in 
this regard. Accuracy is defined according to the relative residual in 
equation ( 3 . 1 . i) , in which r^ = (Z^ G Z) * Pg + Z^g^ and 

is a vector norm in r’^. The minor iteration (see table 1) continues only 

I I r^ 1 I k 

until the required accuracy is attained. Thus, ■ ' ' ■ , < ’"i defines 

l|2'i‘'l|- 

the termination criteria for the minor iterations. 

When the algorithm is far from the solution the reduced gradient — 

1 Iz^g^l 1 — is large and little work is required to locate a direction 
satisfying the acceptance criteria. Only the basic and super basic 
variables are optimized. If one of these variables hits a bound, the 
constraining variable is transferred into the set of nonbasics As 

1 |Z g II is reduced the acceptance criteria becomes more restrictive and 
the current solution to the direction finding problem lies closer to the 
Newton direction. 

At this point, the nonbasic variables must be tested for optimality. 
First order estimates for the Lagrange multipliers are computed as follows: 



n 


-t 


.-1 


lij - K 






s 
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In this environment non-basic variables that reduce the objective function 
when moving away from their bounds (i.e., if < 0 and 

n n 

> 0 and along with free non-basic variables (i.e., u^) are 

n nn n,n“ 

called eligible. Eligible variables are transferred to the superbasic set 

[x ] in conjunction with a maximal basis [A], and the TN algorithm continues 

the next niajor iteration with the new partition. 

A sizable portion of the algorithm’s execution time involves computing 

the search direction p^. While the truncated-Newton method can use any 

iterative method fol solving equation (3.1.i), we have chosen the linear 

conjugate-gradient [CG] method. Although the reduced Hessian matrix 

Z^GZ is typically dense, the product required by [CG], (Z^GZ)'v, is easily 

computable due to the sparsity of the large-scale components. 

The success of the conjugate-gradient method depends upon locating a 

"good" search direction in a small number of iterations. Thus, preconditioning 

the reduced Hessian by the matrix P is important so as to reduce the number of 

CG iterations. Whereas the usual initial element of the CG sequence is 
-k -k 

g , the vector P * g becomes an initial element when preconditioning, 
s s 

where P is a positive-definite matrix. See [2] for further details. 

3.2 Simplicial Decomposition Algorithm 

While the TN algorithm is capable of solving large nonlinear network 
problems, we felt that a first order approximation would be better suited 
for ultra-large examples — problems with more than 10,000 nodes and 
100,000 arcs. The simplicial decomposition algorithm was selected by us to 
meet this goal. This algorithm is best examined in the context of general 
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decomposition methods. .Following Geoff rion [6] we may place the algorithm, 
in the class of "inner linearizaticn/restriction" methods. The SD algorithm 
iterates between (1) solving a linearized subproblem and (2) solving a 
nonlinear master problem on a restricted space subject to non-negativity 
constraints. Table 2 depicts the overall flow. The algorithm was first 
presented in this form by von Hohenbalken [12,20]. Holloway [8] proposed 
the same algorithm, as an extension of the Frank-Wolfe method, and recently 
Lauphongpanich and Hearn [10] devised a restricted verion for the traffic 
assignment problem. We have developed a version of SD specialized to 
handle [NLGN], whereby the master problem is solved Inexactly [13]. 

The following theorem due to Caratheodory provides the necessary 
theoretical foundation for the algorithm. 

Theorem : 

Let X e r’^ be a non-empty convex polytope. Then every 
X e X lies in the relative interior of one of a finite 
number of simplices whose vertices are extreme points of X. 

See [19] for a proof. The main idea behind SD is simple in principle: 

(1) First solve a linearized subproblem to get the extreme points of X. 

We need not generate all the extreme points of the feasible region; 
this would result in a problem as complex as [NLGN]. Instead we 
generate extreme points as needed — in a manner reminiscent of 
Dantzing's column generation method. A (K-1) dimensional simplex 
is defined from K extreme points, and the search for the optimum 
is now restricted on the generated simplex. 




Table 2 ; The SinpHcial Decomposhioo Algorithm 









(2) The master problem attempts to optimize the original function on the 

simplex generated by the subproblem. We have 

Minimize F(w y) 
w 

Subject to: 

K 

J w = 1 

i=l ^ 

0 < w . < 1 
— 1 — 

where y = ,y 2 . • • • tyj^ is the set of generated extreme points and 

w = Wj^,W 2 « • • • sre associated weights. We may reduce further the 

dimension of the master problem making use of the implicit function 
theorem: 

K-1 

Minimize F( w^(y^ - yj^)) 
i=l 


Subject to: 

0 < w. i = 1 , . . . ,K-1 

— 1 . 

Thus we have a nonlinear problem in dimension (K-1) subject to simple non- 
negativity constraints. 

Reducing the master problem to a sequence of unrestricted (K-1) 
dimensional problems, this problem can be solved by any unconstrained algorithm. 
Refer to the work by Mulvey et al. [13] for computational experiments. The 
solution to the master problem is more important when the solution to [NLGN] 
lies in the current simplex. Thus we adjust in a dynamic fashion the 
accuracy in solving the master problem. Again a forcing sequence 
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> 

{r, } — > 0 is employed and the master problem terminates when 
1 * sl I ^ T", — where D is the derived linear basis representing the 

current simplex. Once this degree of accuracy is achieved we return to 
the subproblem. 

Simplicial decomposition provides us with a modular algorithm for 
solving [NLGN]. It converges rapidly to a good approximate solution, 
represented as a linear combination of a few extreme points (K « n) . 

If high accuracy is requires — i.e., exact representation of the optimum 
solution — then we need generate n + 1 extreme points and the master 
problem becomes as difficult as the original problem. In addition the 
subproblem preserves any special structure that may be present in the 
original problem. For [OPT] the subproblem Iterations consists of solving 
a linear generalized network problem. Code LPNETG [12] was employed, 
modified to allow restarting from the basis of the last subproblem. The 
master problem was solved using a Quasi-Newton algorithm. Reference [13] 
provides further details for the implementation of SD and accompanying 
results from computational experiments. 
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4. Data Requirements • 

One of the most important aspects of modeling real vorld systems is 
the ability to collect the required data in a timely fashion. In the case 
of the Models [OPT], data may be classified in two categories: static and 
dynamic. By static we mean data that do not change over a long period of 
time (e.g., airport locations) and by dynamic we refer to data that change 
with time (e.g., flights scheduled), or with technological innovation 
(e.g., aircraft fuel burn rates, navigation systems and flight management). 
For the model to be useful the input data must be readily available. A key 
component of the comprehensive National Airspace System Plan is a centralized 
data-base of aircraft scheduled for, or actually flying the high altitude 
jet routes [7]. This data-base can serve on a real-time basis for updating 
the data required for [OPT]. 

For the prototype model developed the following data were required: 

(I) Airports information 

(II) Flights information 

(III) Fuel burn data. 

Table 3 provides more details. Some data, like the airport coordinates, were 
available through sources used in the past by FAA, \diile other data had to 
be collected for the network model. The following sources were employed: 

(1) International Official Airlines Guide (lOAG) tape, 
providing information about the airports 

(II) Flight Progress Strip data collected by the Control 
Center, providing flights information 

(III) Fuel Burn Model developed by the FAA providing data 

about fuel burn rates for different types of aircrafts. 
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* Airports InformaliOD 

1. Airport ID code 

2. Geographical coordinates 


* Flight Information 

1. Flight ID 

2. Origin airport 

3. Destination airport 

4. Cruise altitude on entering target sector 

5. Cruise altitude on exit from the target sector 

6. Time flight enters the target sector 

7. Time flight exits from the target sector 

8. Flight Hemi code defliung legitimate cruise altitudes 


* Fuel Burn Data 

1. Aircraft type 

2. Fuel burn rate per hour for every legitimate cruise altitude 

3. Fuel burn rate per nautical mile for every legitimate cruise altitude 


Table 3 : Model data requirements 




5. Modeling a Sector of the Indianapolis Control Center 

A network model was built for the airspace centre] lee by a sector of 
the Indianapolis center. The purpose of the model Is to serve as a 
prototype to illustrate the use of the optimization algorithms described 
earlier as well as the feasibility of the proposed model. Data were 
collected for a high traffic period on January 9, 1985 in which a total of 
185 aircrafts crossed the sector over a 6-hour period. The duration of 
flight through the sector ranged from A to 23 minutes. Five distinct cruise 
altitudes above 29,000 feet (FL 290) were selected by the planes. 

The model was built as a multi-period network as shown in figure 3. 

The following provisions were made in the model: 

(I) Allow for delays at the origin airport, up' to three 10-minute 
intervals and similar delays at the destination airports. 

This time grid can be made finer by considering a larger 
number of progressively smaller delay intervals (e.g., six 
5-minute intervals). The added accuracy will be balanced by 
the larger network that has to be solved. 

(II) Allow for every plane to follow one or two alternative cruise 
altitudes besides the one currently followed. Choice was 
restricted to the cruise altitudes one level above and one 
level below the primary altitude. Again, this restriction can 
be relaxed at the expense of generating larger network models — 
the aircrafts could be instructed to follow any one of the four 
or five legitimate cruise altitudes, specified for the 
particular flight. 




FiRorc 3 : rmlofype Nftw«r% Mwlfl for M«n«Rrmrol 
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The model includes the dual objective of assessing risk and cost, 
as proposed in Section 2. Delay cost was considered as a function of 
fuel burn data. Other delay costs like crew salaries and a surrogate 
for customer dissatisfaction could also be included. The risk analysis 
was based on occupancy rates (congestion) at the same altitude during 
time intervals of interest. 

We define the occupancy at level L during a period T as: 

Time spent by aircraft on route i on the same 
LT ^ altitude (L) as aircraft on route j during period T 
ij Time spent by aircraft on route i 

altitude L during period T 


The system-wide relative risk was then defined as: 


I I I 

T L i.j 





■X.) 


2 


where x^,x^ indicate the number of flights on arcs i and j respectively. 
Summation was taken over all legitimate cruise altitudes (L) , and over ten 
36-minute intervals that cover the six hour planning period. Higher 
Interactions — between more than two planes — were ignored. Figure 4 
depicts the different phases of the modeling procedure. 

Example ; 

(1) If a single plane travels at altitude 39000 feet, and the plane 

is in the target section during the whole interval of T = 30 min., 
its contribution to the risk function is 0. 

(2) If a second aircraft flies at 30000 feet for T = 30 min., the 
corresponding risk coefficient is 1. 

(3) If the second plane stays in the target sector for 15 min., 
the risk coefficient is 0.5. 
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Figiire 4; Modeling the Indianapolis Control Sector 
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This formulation provides for two degrees of freedom in the decision 
making process: delay the departure/arrival times in an optimal fashion 
and instruct aircraft to follow alternative altitudes to reduce congestion. 

The resulting network problem has 1295 nodes, 2873 arcs and approximately 
15000 non-zero coefficients describing interacting flights. The test case 
was solved using code NLPNETG. By varying relative weights on the 
transportation and risk function in a systematic fashion the risk/cost 
efficient frontier was traced (figure 5). Again, the efficient frontier 
is not meant to serve as a direct way of comparing risk with cost. Instead 
it guides one in evaluating alternative modes of operation of the air- 
traffic control system, as generated by the model, or with currently 
followed procedures. 

The major advantage of this methodology is that it generates a 
sequence of alternatives that are efficient; i.e., both risk and cost 
values cannot be improved simultaneously. This is easier to understand 
if we notice the location of point A in figure 5 — this point was obtained 

by solving the optimization problem inexactly. From point A we may move to 

( 

a series of alternative solutions for which the system is better off, both 
with respect to transportation cost and risk. 

To study the effect of airplane congestion, we developed a histogram of 
all planes flying at a particular altitude, during the ten time intervals of 
interest. Figure 6 summarizes the results for three particular altitudes, 
before and after the optimization model was used. Note that as expected planes 
were diverted from a highly congested altitude (35000 ft.) to less congested 
routes (31000 ft. and 39000 ft.). This result was obtained with relative 
weights 0.5/9. 5 on both risk and transportation costs. 




-3 
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6. Computational Results 

The FAA Control model can grow in size to challenge the capabilities 
of general purpose algorithms and conventional mainframe computer systems. 

The developed Indianapolis prototype consists of 1295 nodes and 2873 arcs. 

If instead we had chosen six 10-minute intervals for possible delays, and 
five alternative cruise altitudes the resulting problem would consist of 
2A05 nodes and 8510 arcs. If we consider simultaneously ten control sectors, 
with the same number of flights as the Indianapolis center, the network 
grows to approximately 20000 nodes and 80000 arcs. 

To demonstrate the efficiency of the algorithms described in Section 3, 
we have solved several problems arising from a wide range of applications. 
Table A presents relevant information concerning the test problems. First, 
eighteen test problems were solved with the general purpose code MINOS llA] . 
The results from this program formed a benchmark against which to compare 
the truncated Newton — code NLPNETG — and the simplicial decomposition — 
code NGSD — algorithms. Problems that cannot be solved efficiently with 
the general purpose code can be solved with minimal computational resources 
using the specialized network algorithms. Tables 5 and 6 summarize the 
results. We observe that efficiency of the specialized network codes 
Increases with problem size. 

Advances in parallel processing computers are expected to ensure the 
feasibility of new applications. Specifically the air-traffic control 
model will benefit from the use of supercomputers in two domains: 
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(i) Efficiently solve larger models that cover more than one 
control center and ultimately extent to cover the whole 
continental U.S. airspace, with finer time discretization 
(ii) Solve the operational planning model, under real-time 
conditions. 

To illustrate the situation we have’ specialized NLPNETG for the CRAY-XMP 
vector computer at Boeing Computer Services. The optimizing compiler 
proved to be only marginally effective, due to the sparsity of the net- 
work problems. We had to analyze the algorithm in a way to take advantage 
of the pipeline features, and the presence of multiple vector functional 
units. Table 7 summarizes some of the results. We observe that much 
can be gained by specializing the network algorithm for the architecture 
of vector supercomputers. Comparisons with a VAX 11/750 and an IBM 3081 
large mainframe are highlighted in table 8. 

Finally the network model was solved for a range of relative weights, 
with code NLPNETG. We observe (table 9) that a complete analysis can be 
carried out within a few hours, even on a VAX minicomputer. Giving more 
emphasis to the risk function (nonlinear component) causes the problem to 
become, algorithmically, more difficult. This difficulty reflects the 
price we have to pay in going from a linear, transportation cost minimizing 
model, to a nonlinear model that incorporates the nonlinear form of risk 


minimization. 
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7. Conclusions and Future Research 

We have discussed in this paper a general optimization planning model 
for the U.S. air-traffic system. The use of this model has been demonstrated 
using air-traffic data from a representative sector. By developing 
specialized algorithms, we were able to solve the resulting problem 
efficiently, thus demonstrating the feasibility of the approach for strategic 
planning. In addition, taking advantage of the latest technological develop- 
ments in supercomputer design we were, able to solve very large probleems in a 
matter of seconds; thus, the network model can be used as an operational 
planning tool under real-time conditions. 

This research has established that network models can be used as the 
basis for assessing some aspects of the U.S. airspace. The technology — 
in terms of computer systems and algorithms — is available, and the required 
data can be collected. Further research is needed, however, in modeling 
the air-traffic system. We have considered a deterministic model of risk. 

In general a stochastic analysis which also considers the systemwide optimiza- 
tion aspects would be more appropriate. While congestion was used as a risk 
surrogate, other criteria like expected number of aircraft conflicts, or 
expected number of controller intervention should be evaluated as alternatives 
This is an area of potential interface between the optimization model 
described here and the probabilistic models of Odoni and Endoh [15] . 

Other aspects of the air-traffic control system can be examined within 
a network optimization framework, such as the flow control problem [1] , or 
aircraft scheduling problems arising from aggregate solutions of the current 
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model. Issues like the uncertainty involved in determining projected 
demand for air transportation merit investigation, and the end effects 
due to the finite planning horizon have to be examined. As another 
extension, resource limitations can be Imposed by Introducing arc 
multipliers on traffic arriving/departing from an airport, thus controlling 
the total number of passengers an airport facility can handle. 

While substantial progress has been made towards building and 
verifying the model, additional work needs to be done in model validation. 
Alternative strategies generated by the model have to be compared with 
existing methodologies to establish the correspondence of the model and 
its results to the perceived reality. This is another area of potential 
interface between the optimization model and the general probabilistic 
framework presented in Powell, Mulvey and Babu [16]. 
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^IGINAL PAGE fS 
W POOR DUALITY 


PROBLEM 

Size 

(Nodes/Arcs) 

Free arcs 
at optimum 

Condition No. 
of Reduced 
Hessian 

Objective 

value 

I« Norm 
of Reduced 
Gradient 

Description | 

PTN30 

30/ 46 

15 

>10* 


.0045 

Dallas Water 

PTN150 

150/ 196 

44 

>10* 

..4819730E5 . 

.0020 

Distribution 

PTN660 

666/ 906 

240 

>10® 

-.2061074E6 

.0160 

models 

SMBANK 

64/ 117 

54 

>10* 

•.7129290E7 

.0003 

Matrix balanc- 
ing models 

1 

BIGBANK 

1116/2230 

946 

>10* 

-.4205693E7 

.0004 

RANBANK 

800/2000 

13 

>10* 

.3018788E6 

.4025 

STICK 1 

209/ 454 

246 

^KP 

.6934392E1 

.0001 

Stick percula- i 

STICK2 

650/1412 

763 

>10* 

.3124563E1 

.0001 

tion models, { 

STICKS 

782/1686 

905 

^10* 

.1117978E2 

.0006 

electrical net- ' 

ST1CK4 

832/2264 

1433 

>10* 

.I56619SE1 

.0010 

works ! 

GROUPJac 

200/ 500 

100 

>10* 

.1011792E5 

.0060 

Randomly gen- 

GROUPlad 

400/1000 

119 

>10* 

.3834884E5 

.0049 

crated, strictly 

GROUPlae 

800/2000 

230 

>10* 

.22657 14E6 

.4070 

convex net- 

works 

MARK] 

23/ 42 

4 

>10* 

-.1145214E6 


Markowitz 

MARK2 

53/ 102 

10 

>10* 

-.3271077E6 

.2420 

protfolio con- 

MARK3 

857/1710 

35 

>10* 

.1150057E8 

.0710 

struction 

FAATEST 

28/ 64 

6 

>10* 

.7725644E5 


Air traffic con- 

FAA.1ND3 

1295/2873 

787 

>10* 

.2S56325E1 


trol model 


Table 4 : Te*t Problems 


Problem 

NLPNETG Solution 

MINOS Solution 

Time (sec) 

KiKCSS3l 

Time (sec) 

BRSZBII 

PTN30 

7.33 

4.5E-3 

44.08 

8.8E-6 

PIN 1 50 

23.86 

2.0E-3 

305.01 

8.8E-6 

PTN660 

297.93 

I.6E-2 

9647.35 

4.0E-4 

SMBANK 

21.50 

3.0E-4 

287.20 

4.9E-4 

BIGBANK 

9100.00 

4.0E-4 

28800.00* 

4.0E-1 

RANBANK 

245.37 

4.0E-1 

5477.97* 

l.OE+2 

STICK 1 

19.05 

l.OE-4 

12392.00 

3.9E-6 

STJCK2 

103.12 

J.OE-4 

»• 


STICKS 

73.75 

6.0E-4 

•• 


ST1CK4 

166.50 

l.OE-3 

•• 


GROUPlae 

1652.00 

6.0E-3 

2312.69 

5.9E-2 

GROUPlad 

10227.00 

4.9E-3 

20347.00 

5.9E-1 

GROUPlae 

7376.00 

4.1E-1 

16115.00 

2.8E-1 

MARKl 

3.07 

2.0E-3 

12.23 

1.3E-1 

MARK2 

21.60 

2.4E-1 

33.03 

1.4E-1 

MARK 3 

204.43 

7.1 E-2 

1341.53 

1.5E-1 

FAA.TEST 

U4 

5.0E-2 

15.22 

3.4E-9 

FAA.IND3 

707.00 

5.0E-2 

111600.00 

2.8E-4 


* Ttmixui^d wbcB 0J03%fvocBPptUMB> 


** Did bot cetmiyp ifier I boun 


Table 5 : Comparison of NLPNETG with MINOS 
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Problem 

NGSD Solution 

MINOS Solution 

Time (sec) 


Time (sec) 

BRSSQi 

PTN30 



44.08 

8.8E-6 

PTN150 



305.01 

8.8E-6 

PTN660 



9647.35 

4.0E-4 

RANBANK 


4.0E-1 

5477.97* 

l.OE+2 

STICK 1 

50.25 

0.40 

12392.00 

3.9E-6 

ST1CK2 

517.48 

0.94 



STICK3 

1387.25 

0.08 

00 


ST1CK4 

2688.72 

0.15 

00 


GROUPlac 

314.85 

1.60 

2312.69 

5.9E-2 

GROUPlad 

» 606.90 

0.40 

20347.00 

5.9E-1 

GROUPlac 

395.38 

0.52 

16115.00 

2.8E-1 

MARKl 

3.07 

2.0E-3 

12.23 

1.3E-1 

MARK2 

21.60 

2.4E-1 

33.03 

1.4E-1 

MARK 3 

204.43 

7.1E-2 

1341.53 

1.5E-1 


* T«nni7^Aird when within 0X)2% from Optimum 
** Did Dot converge nftci 6 boun 


Tabic 6 : Comparison of GNSD with MINOS 


Problem 

NLPNETG Solution times (sec) 


Without vectorization 

Compiler vectorization 

User vectorization 

PTN 1 50 

0.328 


0.165 

PTN660 

2.435 


1.402 

SMBANK 

0.358 


0.177 

BIGBANK 

244.107 

225.900 

58.896 

STICK 4 

0 

0 

2.925 

GROUPlac 

18.668 

17.215 

4.983 

GROUPlac 

0 

• 

49.454 

MARK 3 

• 

0 

2.131 


* Problem oot rolwed with thie opboo 


Tabic 7 : Vectorization of NLPNETG on the CRAY/XMP 
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pace is 
OUALITV 


Problem 

NL 

PNETG Solution times (sec) | 

IBM 308 1 

VAX 11/750 (Unix) 

CRAY/XMP 

PTN150 

1.37 

23.86 


PTN660 

36.68 

297.93 

1.402 

SMBANK 

2,33 

21.50 

0.177 

BIGBANK 

884.53 

9100.00 

58.896 

GROUPlac 

218,46 

1652,00 

4.983 

GROUP Jac 

1320.82 

10227.00 

49.454 

MARK3 

24.87 

204.43 

2.131 

Average 

281.19(17) 

3075.25 (184) 

16.744 (1) 


Table 8 : Testing NLPNETG on difTcrent computer systems 


Weight ( H j ) 

Time (sec) 

/„ norm 

J.O 

J2JJ 

0.020 

0.999999 


0.054 

0.999997 

396 

0.094 

0.999995 

396 

0.150 

0;999992 

406 

0.187 

0.999990 

367 

0.467 

0.999975 

357 

0.069 

0.0 

18 

0.0 


Table 9 : Tracing the efHcient frontier with NLPNETG 
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